home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 2002 November / SGI IRIX Base Documentation 2002 November.iso / usr / share / catman / p_man / cat3 / SCSL / dlaebz.z / dlaebz
Encoding:
Text File  |  2002-10-03  |  12.9 KB  |  331 lines
  1.  
  2.  
  3.  
  4. DDDDLLLLAAAAEEEEBBBBZZZZ((((3333SSSS))))                                                          DDDDLLLLAAAAEEEEBBBBZZZZ((((3333SSSS))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      DLAEBZ - contain the iteration loops which compute and use the function
  10.      N(w), which is the count of eigenvalues of a symmetric tridiagonal matrix
  11.      T less than or equal to its argument w
  12.  
  13. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  14.      SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, RELTOL,
  15.                         PIVMIN, D, E, E2, NVAL, AB, C, MOUT, NAB, WORK, IWORK,
  16.                         INFO )
  17.  
  18.          INTEGER        IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
  19.  
  20.          DOUBLE         PRECISION ABSTOL, PIVMIN, RELTOL
  21.  
  22.          INTEGER        IWORK( * ), NAB( MMAX, * ), NVAL( * )
  23.  
  24.          DOUBLE         PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( *
  25.                         ), WORK( * )
  26.  
  27. IIIIMMMMPPPPLLLLEEEEMMMMEEEENNNNTTTTAAAATTTTIIIIOOOONNNN
  28.      These routines are part of the SCSL Scientific Library and can be loaded
  29.      using either the -lscs or the -lscs_mp option.  The -lscs_mp option
  30.      directs the linker to use the multi-processor version of the library.
  31.  
  32.      When linking to SCSL with -lscs or -lscs_mp, the default integer size is
  33.      4 bytes (32 bits). Another version of SCSL is available in which integers
  34.      are 8 bytes (64 bits).  This version allows the user access to larger
  35.      memory sizes and helps when porting legacy Cray codes.  It can be loaded
  36.      by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
  37.      only one of the two versions; 4-byte integer and 8-byte integer library
  38.      calls cannot be mixed.
  39.  
  40. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  41.      DLAEBZ contains the iteration loops which compute and use the function
  42.      N(w), which is the count of eigenvalues of a symmetric tridiagonal matrix
  43.      T less than or equal to its argument w. It performs a choice of two types
  44.      of loops:
  45.  
  46.      IJOB=1, followed by
  47.      IJOB=2: It takes as input a list of intervals and returns a list of
  48.              sufficiently small intervals whose union contains the same
  49.              eigenvalues as the union of the original intervals.
  50.              The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
  51.              The output interval (AB(j,1),AB(j,2)] will contain
  52.              eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
  53.  
  54.      IJOB=3: It performs a binary search in each input interval
  55.              (AB(j,1),AB(j,2)] for a point  w(j)  such that
  56.              N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
  57.              the search.  If such a w(j) is found, then on output
  58.              AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
  59.              (AB(j,1),AB(j,2)] will be a small interval containing the
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. DDDDLLLLAAAAEEEEBBBBZZZZ((((3333SSSS))))                                                          DDDDLLLLAAAAEEEEBBBBZZZZ((((3333SSSS))))
  71.  
  72.  
  73.  
  74.              point where N(w) jumps through NVAL(j), unless that point
  75.              lies outside the initial interval.
  76.  
  77.      Note that the intervals are in all cases half-open intervals, i.e., of
  78.      the form  (a,b] , which includes  b  but not  a .
  79.  
  80.      To avoid underflow, the matrix should be scaled so that its largest
  81.      element is no greater than  overflow**(1/2) * underflow**(1/4) in
  82.      absolute value.  To assure the most accurate computation of small
  83.      eigenvalues, the matrix should be scaled to be
  84.      not much smaller than that, either.
  85.  
  86.      See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal Matrix",
  87.      Report CS41, Computer Science Dept., Stanford
  88.      University, July 21, 1966
  89.  
  90.      Note: the arguments are, in general, *not* checked for unreasonable
  91.      values.
  92.  
  93.  
  94. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  95.      IJOB    (input) INTEGER
  96.              Specifies what is to be done:
  97.              = 1:  Compute NAB for the initial intervals.
  98.              = 2:  Perform bisection iteration to find eigenvalues of T.
  99.              = 3:  Perform bisection iteration to invert N(w), i.e., to find a
  100.              point which has a specified number of eigenvalues of T to its
  101.              left.  Other values will cause DLAEBZ to return with INFO=-1.
  102.  
  103.      NITMAX  (input) INTEGER
  104.              The maximum number of "levels" of bisection to be performed,
  105.              i.e., an interval of width W will not be made smaller than 2^(-
  106.              NITMAX) * W.  If not all intervals have converged after NITMAX
  107.              iterations, then INFO is set to the number of non-converged
  108.              intervals.
  109.  
  110.      N       (input) INTEGER
  111.              The dimension n of the tridiagonal matrix T.  It must be at least
  112.              1.
  113.  
  114.      MMAX    (input) INTEGER
  115.              The maximum number of intervals.  If more than MMAX intervals are
  116.              generated, then DLAEBZ will quit with INFO=MMAX+1.
  117.  
  118.      MINP    (input) INTEGER
  119.              The initial number of intervals.  It may not be greater than
  120.              MMAX.
  121.  
  122.      NBMIN   (input) INTEGER
  123.              The smallest number of intervals that should be processed using a
  124.              vector loop.  If zero, then only the scalar loop will be used.
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. DDDDLLLLAAAAEEEEBBBBZZZZ((((3333SSSS))))                                                          DDDDLLLLAAAAEEEEBBBBZZZZ((((3333SSSS))))
  137.  
  138.  
  139.  
  140.      ABSTOL  (input) DOUBLE PRECISION
  141.              The minimum (absolute) width of an interval.  When an interval is
  142.              narrower than ABSTOL, or than RELTOL times the larger (in
  143.              magnitude) endpoint, then it is considered to be sufficiently
  144.              small, i.e., converged.  This must be at least zero.
  145.  
  146.      RELTOL  (input) DOUBLE PRECISION
  147.              The minimum relative width of an interval.  When an interval is
  148.              narrower than ABSTOL, or than RELTOL times the larger (in
  149.              magnitude) endpoint, then it is considered to be sufficiently
  150.              small, i.e., converged.  Note: this should always be at least
  151.              radix*machine epsilon.
  152.  
  153.      PIVMIN  (input) DOUBLE PRECISION
  154.              The minimum absolute value of a "pivot" in the Sturm sequence
  155.              loop.  This *must* be at least  max |e(j)**2| * safe_min  and at
  156.              least safe_min, where safe_min is at least the smallest number
  157.              that can divide one without overflow.
  158.  
  159.      D       (input) DOUBLE PRECISION array, dimension (N)
  160.              The diagonal elements of the tridiagonal matrix T.
  161.  
  162.      E       (input) DOUBLE PRECISION array, dimension (N)
  163.              The offdiagonal elements of the tridiagonal matrix T in positions
  164.              1 through N-1.  E(N) is arbitrary.
  165.  
  166.      E2      (input) DOUBLE PRECISION array, dimension (N)
  167.              The squares of the offdiagonal elements of the tridiagonal matrix
  168.              T.  E2(N) is ignored.
  169.  
  170.      NVAL    (input/output) INTEGER array, dimension (MINP)
  171.              If IJOB=1 or 2, not referenced.  If IJOB=3, the desired values of
  172.              N(w).  The elements of NVAL will be reordered to correspond with
  173.              the intervals in AB.  Thus, NVAL(j) on output will not, in
  174.              general be the same as NVAL(j) on input, but it will correspond
  175.              with the interval (AB(j,1),AB(j,2)] on output.
  176.  
  177.      AB      (input/output) DOUBLE PRECISION array, dimension (MMAX,2)
  178.              The endpoints of the intervals.  AB(j,1) is  a(j), the left
  179.              endpoint of the j-th interval, and AB(j,2) is b(j), the right
  180.              endpoint of the j-th interval.  The input intervals will, in
  181.              general, be modified, split, and reordered by the calculation.
  182.  
  183.      C       (input/output) DOUBLE PRECISION array, dimension (MMAX)
  184.              If IJOB=1, ignored.  If IJOB=2, workspace.  If IJOB=3, then on
  185.              input C(j) should be initialized to the first search point in the
  186.              binary search.
  187.  
  188.      MOUT    (output) INTEGER
  189.              If IJOB=1, the number of eigenvalues in the intervals.  If IJOB=2
  190.              or 3, the number of intervals output.  If IJOB=3, MOUT will equal
  191.              MINP.
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. DDDDLLLLAAAAEEEEBBBBZZZZ((((3333SSSS))))                                                          DDDDLLLLAAAAEEEEBBBBZZZZ((((3333SSSS))))
  203.  
  204.  
  205.  
  206.      NAB     (input/output) INTEGER array, dimension (MMAX,2)
  207.              If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).  If
  208.              IJOB=2, then on input, NAB(i,j) should be set.  It must satisfy
  209.              the condition:  N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
  210.              which means that in interval i only eigenvalues
  211.              NAB(i,1)+1,...,NAB(i,2) will be considered.  Usually,
  212.              NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with IJOB=1.
  213.              On output, NAB(i,j) will contain
  214.              max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of the
  215.              input interval that the output interval (AB(j,1),AB(j,2)] came
  216.              from, and na(k) and nb(k) are the the input values of NAB(k,1)
  217.              and NAB(k,2).  If IJOB=3, then on output, NAB(i,j) contains
  218.              N(AB(i,j)), unless N(w) > NVAL(i) for all search points  w , in
  219.              which case NAB(i,1) will not be modified, i.e., the output value
  220.              will be the same as the input value (modulo reorderings -- see
  221.              NVAL and AB), or unless N(w) < NVAL(i) for all search points  w ,
  222.              in which case NAB(i,2) will not be modified.  Normally, NAB
  223.              should be set to some distinctive value(s) before DLAEBZ is
  224.              called.
  225.  
  226.      WORK    (workspace) DOUBLE PRECISION array, dimension (MMAX)
  227.              Workspace.
  228.  
  229.      IWORK   (workspace) INTEGER array, dimension (MMAX)
  230.              Workspace.
  231.  
  232.      INFO    (output) INTEGER
  233.              = 0:       All intervals converged.
  234.              = 1--MMAX: The last INFO intervals did not converge.
  235.              = MMAX+1:  More than MMAX intervals were generated.
  236.  
  237. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  238.          This routine is intended to be called only by other LAPACK routines,
  239.      thus the interface is less user-friendly.  It is intended for two
  240.      purposes:
  241.  
  242.      (a) finding eigenvalues.  In this case, DLAEBZ should have one or
  243.          more initial intervals set up in AB, and DLAEBZ should be called
  244.          with IJOB=1.  This sets up NAB, and also counts the eigenvalues.
  245.          Intervals with no eigenvalues would usually be thrown out at
  246.          this point.  Also, if not all the eigenvalues in an interval i
  247.          are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
  248.          For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
  249.          eigenvalue.  DLAEBZ is then called with IJOB=2 and MMAX
  250.          no smaller than the value of MOUT returned by the call with
  251.          IJOB=1.  After this (IJOB=2) call, eigenvalues NAB(i,1)+1
  252.          through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
  253.          tolerance specified by ABSTOL and RELTOL.
  254.  
  255.      (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
  256.          In this case, start with a Gershgorin interval  (a,b).  Set up
  257.          AB to contain 2 search intervals, both initially (a,b).  One
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268. DDDDLLLLAAAAEEEEBBBBZZZZ((((3333SSSS))))                                                          DDDDLLLLAAAAEEEEBBBBZZZZ((((3333SSSS))))
  269.  
  270.  
  271.  
  272.          NVAL element should contain  f-1  and the other should contain  l
  273.          , while C should contain a and b, resp.  NAB(i,1) should be -1
  274.          and NAB(i,2) should be N+1, to flag an error if the desired
  275.          interval does not lie in (a,b).  DLAEBZ is then called with
  276.          IJOB=3.  On exit, if w(f-1) < w(f), then one of the intervals --
  277.          j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
  278.          if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
  279.          >= 0, then the interval will have  N(AB(j,1))=NAB(j,1)=f-k and
  280.          N(AB(j,2))=NAB(j,2)=f+r.  The cases w(l) < w(l+1) and
  281.          w(l-r)=...=w(l+k) are handled similarly.
  282.  
  283.  
  284. SSSSEEEEEEEE AAAALLLLSSSSOOOO
  285.      INTRO_LAPACK(3S), INTRO_SCSL(3S)
  286.  
  287.      This man page is available only online.
  288.  
  289.  
  290.  
  291.  
  292.  
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299.  
  300.  
  301.  
  302.  
  303.  
  304.  
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.                                                                         PPPPaaaaggggeeee 5555
  328.  
  329.  
  330.  
  331.